Col Ext by Region

library(data.table)
library(ggplot2)
library(lme4)
library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:data.table’:

    between, first, last

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

load("~/Documents/grad school/Rutgers/LabWork/Colonization_Extinction/spp_master_ztemp.Rdata")

First, let’s get some regional statistics

#count colonizations by region
col <- spp_master_ztemp %>%
  group_by(reg) %>%
  summarise(sum(col))
#count extinctions by region
ext <- spp_master_ztemp %>%
  group_by(reg) %>%
  summarise(sum(now_ext))
#count observations by region
count <- spp_master_ztemp %>%
  group_by(reg) %>%
  summarise(n())
#first year
year.first <- spp_master_ztemp %>%
  group_by(reg) %>%
  summarise(min(year))
#last year
year.last <- spp_master_ztemp %>%
  group_by(reg) %>%
  summarise(max(year))
region_summaries <- cbind(count, col[,2], ext[,2], year.first[,2], year.last[,2])
colnames(region_summaries) <- c("reg", "n", "n.col", "n.ext", "year.first","year.last")

Now we’ll visualize regions in a plot

ggsave(plot = g, filename = "plots/occurence_plot.png")

Saving 7 x 7 in image

g

Now, let’s do the same process as we did in 04_24_19_Temp_link_traits_fish BUT region specific

Surface colonization comparison by region

Now, I’ll merge all surface colonization models

surface_colonization_model_comparison_byregion <- rbind(surface_colonization_model_comparison_ai, surface_colonization_model_comparison_ebs, surface_colonization_model_comparison_gmex, surface_colonization_model_comparison_goa, surface_colonization_model_comparison_neus, surface_colonization_model_comparison_newf, surface_colonization_model_comparison_sa, surface_colonization_model_comparison_shelf, surface_colonization_model_comparison_wctri)

Best surface colonization models per region

surface_colonization_model_comparison_byregion %>%
  filter(converge == TRUE, !grepl("mean", variable)) %>% #ignoring mean for now, and get rid of those that didn't converge
  group_by(reg) %>% #group by region
  arrange(AICc) %>% #arrange by AICc
  top_n(-5, AICc) #take 5 best models (lowest AICc)

NA

---

Bottom colonization comparison by region

  for (i in 1:length(bottom_variables)){
    mod <- glmer(col ~ get(bottom_variables[i]) + (1|year_factor), family = binomial, data = spp_master_ztemp[reg == "wctri"])
    bottom_colonization_model_comparison_wctri[i,variable := bottom_variables[i]]
    bottom_colonization_model_comparison_wctri[i,coef := coef(summary(mod))[,"Estimate"][2]]
    bottom_colonization_model_comparison_wctri[i,p_value := coef(summary(mod))[,"Pr(>|z|)"][2]]
    bottom_colonization_model_comparison_wctri[i,AICc := AICc(mod)]
    
    print(paste(i, length(bottom_variables), sep = "/"))
      
  }

[1] "1/124"
[1] "2/124"
[1] "3/124"
[1] "4/124"
[1] "5/124"
[1] "6/124"
[1] "7/124"
[1] "8/124"
[1] "9/124"
[1] "10/124"
[1] "11/124"
[1] "12/124"
[1] "13/124"

Model failed to converge with max|grad| = 0.0308063 (tol = 0.001, component 1)

[1] "14/124"
[1] "15/124"
[1] "16/124"
[1] "17/124"
[1] "18/124"
[1] "19/124"
[1] "20/124"
[1] "21/124"
[1] "22/124"

Model failed to converge with max|grad| = 0.0296878 (tol = 0.001, component 1)

[1] "23/124"
[1] "24/124"
[1] "25/124"
[1] "26/124"
[1] "27/124"
[1] "28/124"
[1] "29/124"
[1] "30/124"
[1] "31/124"
[1] "32/124"

Model failed to converge with max|grad| = 0.0298623 (tol = 0.001, component 1)

[1] "33/124"
[1] "34/124"
[1] "35/124"
[1] "36/124"
[1] "37/124"
[1] "38/124"
[1] "39/124"
[1] "40/124"
[1] "41/124"
[1] "42/124"
[1] "43/124"
[1] "44/124"
[1] "45/124"
[1] "46/124"
[1] "47/124"
[1] "48/124"

Model failed to converge with max|grad| = 0.0288633 (tol = 0.001, component 1)

[1] "49/124"
[1] "50/124"

Model failed to converge with max|grad| = 0.0268189 (tol = 0.001, component 1)

[1] "51/124"

Model failed to converge with max|grad| = 0.0266311 (tol = 0.001, component 1)

[1] "52/124"
[1] "53/124"
[1] "54/124"
[1] "55/124"

Model failed to converge with max|grad| = 0.0268531 (tol = 0.001, component 1)

[1] "56/124"
[1] "57/124"

Model failed to converge with max|grad| = 0.028596 (tol = 0.001, component 1)

[1] "58/124"
[1] "59/124"

Model failed to converge with max|grad| = 0.0282459 (tol = 0.001, component 1)

[1] "60/124"
[1] "61/124"
[1] "62/124"
[1] "63/124"
[1] "64/124"

Model failed to converge with max|grad| = 0.0304929 (tol = 0.001, component 1)

[1] "65/124"
[1] "66/124"
[1] "67/124"
[1] "68/124"
[1] "69/124"
[1] "70/124"
[1] "71/124"
[1] "72/124"
[1] "73/124"
[1] "74/124"
[1] "75/124"

Model failed to converge with max|grad| = 0.0316174 (tol = 0.001, component 1)

[1] "76/124"
[1] "77/124"

Model failed to converge with max|grad| = 0.0261989 (tol = 0.001, component 1)

[1] "78/124"
[1] "79/124"
[1] "80/124"
[1] "81/124"
[1] "82/124"
[1] "83/124"
[1] "84/124"
[1] "85/124"
[1] "86/124"
[1] "87/124"
[1] "88/124"

Model failed to converge with max|grad| = 0.0309174 (tol = 0.001, component 1)

[1] "89/124"
[1] "90/124"
[1] "91/124"
[1] "92/124"

Model failed to converge with max|grad| = 0.029357 (tol = 0.001, component 1)

[1] "93/124"
[1] "94/124"
[1] "95/124"
[1] "96/124"
[1] "97/124"

Model failed to converge with max|grad| = 0.0290942 (tol = 0.001, component 1)

[1] "98/124"
[1] "99/124"
[1] "100/124"
[1] "101/124"
[1] "102/124"
[1] "103/124"

Model failed to converge with max|grad| = 0.0294299 (tol = 0.001, component 1)

[1] "104/124"
[1] "105/124"
[1] "106/124"
[1] "107/124"
[1] "108/124"

Model failed to converge with max|grad| = 0.0313106 (tol = 0.001, component 1)

[1] "109/124"
[1] "110/124"
[1] "111/124"

Model failed to converge with max|grad| = 0.0295465 (tol = 0.001, component 1)

[1] "112/124"

Model failed to converge with max|grad| = 0.0314474 (tol = 0.001, component 1)

[1] "113/124"
[1] "114/124"

Model failed to converge with max|grad| = 0.0247102 (tol = 0.001, component 1)

[1] "115/124"
[1] "116/124"
[1] "117/124"
[1] "118/124"

Model failed to converge with max|grad| = 0.0251905 (tol = 0.001, component 1)

[1] "119/124"
[1] "120/124"
[1] "121/124"
[1] "122/124"
[1] "123/124"
[1] "124/124"

Now, I’ll merge all bottom colonization models

bottom_colonization_model_comparison_byregion <- rbind(bottom_colonization_model_comparison_ai, bottom_colonization_model_comparison_ebs, bottom_colonization_model_comparison_gmex, bottom_colonization_model_comparison_goa, bottom_colonization_model_comparison_neus, bottom_colonization_model_comparison_newf, bottom_colonization_model_comparison_sa, bottom_colonization_model_comparison_shelf, bottom_colonization_model_comparison_wctri)

Best bottom colonization models per region

---

Surface extinction comparison by region

Now, I’ll merge all surface extinction models

surface_extinction_model_comparison_byregion <- rbind(surface_extinction_model_comparison_ai, surface_extinction_model_comparison_ebs, surface_extinction_model_comparison_gmex, surface_extinction_model_comparison_goa, surface_extinction_model_comparison_neus, surface_extinction_model_comparison_newf, surface_extinction_model_comparison_sa, surface_extinction_model_comparison_shelf, surface_extinction_model_comparison_wctri)

Best surface extinction models per region

surface_extinction_model_comparison_byregion %>%
  filter(converge == TRUE, !grepl("mean", variable)) %>% #ignoring mean for now, and get rid of those that didn't converge
  group_by(reg) %>% #group by region
  arrange(AICc) %>% #arrange by AICc
  top_n(-5, AICc) #take 5 best models (lowest AICc)
  

---

Bottom extinction comparison by region

bottom_variables <- colnames(spp_master_ztemp[,169:292])
  
#regions

#ai (1:3, 15, 24, 44, 45, 47, 55, 57 fail to converge)
bottom_extinction_model_comparison_ai <- as.data.table(matrix(nrow = length(bottom_variables)))
bottom_extinction_model_comparison_ai[, variable:=as.factor(V1)][, coef:=as.numeric(V1)][, p_value:=as.numeric(V1)][, AICc:=as.numeric(V1)]
bottom_extinction_model_comparison_ai[, V1 := NULL]
  
  for (i in 1:length(bottom_variables)){
    mod <- glmer(now_ext ~ get(bottom_variables[i]) + (1|year_factor), family = binomial, data = spp_master_ztemp[reg == "ai"])
    bottom_extinction_model_comparison_ai[i,variable := bottom_variables[i]]
    bottom_extinction_model_comparison_ai[i,coef := coef(summary(mod))[,"Estimate"][2]]
    bottom_extinction_model_comparison_ai[i,p_value := coef(summary(mod))[,"Pr(>|z|)"][2]]
    bottom_extinction_model_comparison_ai[i,AICc := AICc(mod)]
    
    print(paste(i, length(bottom_variables), sep = "/"))
      
  }
  #add column for those that didn't converge
  bottom_extinction_model_comparison_ai$converge <- T
  bottom_extinction_model_comparison_ai$converge[c(1:3, 15, 24, 44, 45, 47, 55, 57)] <- F
  
  #add column to designate region 
  bottom_extinction_model_comparison_ai$reg <- "ai"
  
#ebs (23 did not converge)

bottom_extinction_model_comparison_ebs <- as.data.table(matrix(nrow = length(bottom_variables)))
bottom_extinction_model_comparison_ebs[, variable:=as.factor(V1)][, coef:=as.numeric(V1)][, p_value:=as.numeric(V1)][, AICc:=as.numeric(V1)]
bottom_extinction_model_comparison_ebs[, V1 := NULL]
  
  for (i in 1:length(bottom_variables)){
    mod <- glmer(now_ext ~ get(bottom_variables[i]) + (1|year_factor), family = binomial, data = spp_master_ztemp[reg == "ebs"])
    bottom_extinction_model_comparison_ebs[i,variable := bottom_variables[i]]
    bottom_extinction_model_comparison_ebs[i,coef := coef(summary(mod))[,"Estimate"][2]]
    bottom_extinction_model_comparison_ebs[i,p_value := coef(summary(mod))[,"Pr(>|z|)"][2]]
    bottom_extinction_model_comparison_ebs[i,AICc := AICc(mod)]
    
    print(paste(i, length(bottom_variables), sep = "/"))
      
  }
  #add column for those that didn't converge
  bottom_extinction_model_comparison_ebs$converge <- T
  bottom_extinction_model_comparison_ebs$converge[c(23)] <- F
  
  #add column to designate region 
  bottom_extinction_model_comparison_ebs$reg <- "ebs"
  
#gmex (5, 6, 7, 8, 13, 15, 17, 18, 20, 21, 23, 24, 25, 27, 28, 31, 32, 34, 36, 40:44, 46:48, 52, 56, 57, 60, 61, 63, 67, 70, 75, 80, 90, 92, 99, 110, 113,  fail to converge)

bottom_extinction_model_comparison_gmex <- as.data.table(matrix(nrow = length(bottom_variables)))
bottom_extinction_model_comparison_gmex[, variable:=as.factor(V1)][, coef:=as.numeric(V1)][, p_value:=as.numeric(V1)][, AICc:=as.numeric(V1)]
bottom_extinction_model_comparison_gmex[, V1 := NULL]
  
  for (i in 1:length(bottom_variables)){
    mod <- glmer(now_ext ~ get(bottom_variables[i]) + (1|year_factor), family = binomial, data = spp_master_ztemp[reg == "gmex"])
    mod1 <- glmer(now_ext ~ get(bottom_variables[i]) + (1|year_factor), family = binomial, data = spp_master_ztemp[reg == "gmex"])
    bottom_extinction_model_comparison_gmex[i,variable := bottom_variables[i]]
    bottom_extinction_model_comparison_gmex[i,coef := coef(summary(mod))[,"Estimate"][2]]
    bottom_extinction_model_comparison_gmex[i,p_value := coef(summary(mod))[,"Pr(>|z|)"][2]]
    bottom_extinction_model_comparison_gmex[i,AICc := AICc(mod)]
    
    print(paste(i, length(bottom_variables), sep = "/"))
      
  }
  #add column for those that didn't converge
  bottom_extinction_model_comparison_gmex$converge <- T
  bottom_extinction_model_comparison_gmex$converge[c(5, 6, 7, 8, 13, 15, 17, 18, 20, 21, 23, 24, 25, 27, 28, 31, 32, 34, 36, 40:44, 46:48, 52, 56, 57, 60, 61, 63, 67, 70, 75, 80, 90, 92, 99, 110, 113)] <- F
  
  #add column to designate region 
  bottom_extinction_model_comparison_gmex$reg <- "gmex"

#goa (62 doesn't converge)

bottom_extinction_model_comparison_goa <- as.data.table(matrix(nrow = length(bottom_variables)))
bottom_extinction_model_comparison_goa[, variable:=as.factor(V1)][, coef:=as.numeric(V1)][, p_value:=as.numeric(V1)][, AICc:=as.numeric(V1)]
bottom_extinction_model_comparison_goa[, V1 := NULL]
  
  for (i in 1:length(bottom_variables)){
    mod <- glmer(now_ext ~ get(bottom_variables[i]) + (1|year_factor), family = binomial, data = spp_master_ztemp[reg == "goa"])
    bottom_extinction_model_comparison_goa[i,variable := bottom_variables[i]]
    bottom_extinction_model_comparison_goa[i,coef := coef(summary(mod))[,"Estimate"][2]]
    bottom_extinction_model_comparison_goa[i,p_value := coef(summary(mod))[,"Pr(>|z|)"][2]]
    bottom_extinction_model_comparison_goa[i,AICc := AICc(mod)]
    
    print(paste(i, length(bottom_variables), sep = "/"))
      
  }
  #add column for those that didn't converge
  bottom_extinction_model_comparison_goa$converge <- T
  bottom_extinction_model_comparison_goa$converge[c(62)] <- T
  
  #add column to designate region 
  bottom_extinction_model_comparison_goa$reg <- "goa"

#neus (1, 4, 9, 10, 12, 16:19, 22, 24, 27, 28, 29, 32, 33, 36, 37, 39, 40, 43,  failed to converge)

bottom_extinction_model_comparison_neus <- as.data.table(matrix(nrow = length(bottom_variables)))
bottom_extinction_model_comparison_neus[, variable:=as.factor(V1)][, coef:=as.numeric(V1)][, p_value:=as.numeric(V1)][, AICc:=as.numeric(V1)]
bottom_extinction_model_comparison_neus[, V1 := NULL]
  
  for (i in 1:length(bottom_variables)){
    mod <- glmer(now_ext ~ get(bottom_variables[i]) + (1|year_factor), family = binomial, data = spp_master_ztemp[reg == "neus"])
    bottom_extinction_model_comparison_neus[i,variable := bottom_variables[i]]
    bottom_extinction_model_comparison_neus[i,coef := coef(summary(mod))[,"Estimate"][2]]
    bottom_extinction_model_comparison_neus[i,p_value := coef(summary(mod))[,"Pr(>|z|)"][2]]
    bottom_extinction_model_comparison_neus[i,AICc := AICc(mod)]
    
    print(paste(i, length(bottom_variables), sep = "/"))
      
  }
#add column for those that didn't converge
bottom_extinction_model_comparison_neus$converge <- T
bottom_extinction_model_comparison_neus$converge[c(1, 4, 9, 10, 12, 16:19, 22, 24, 27, 28, 29, 32, 33, 36, 37, 39, 40, 43)] <- F

#add column to designate region 
  bottom_extinction_model_comparison_neus$reg <- "neus"

#newf (ALL failed to converge (all singularities) )

bottom_extinction_model_comparison_newf <- as.data.table(matrix(nrow = length(bottom_variables)))
bottom_extinction_model_comparison_newf[, variable:=as.factor(V1)][, coef:=as.numeric(V1)][, p_value:=as.numeric(V1)][, AICc:=as.numeric(V1)]
bottom_extinction_model_comparison_newf[, V1 := NULL]
  
  for (i in 1:length(bottom_variables)){
    mod <- glmer(now_ext ~ get(bottom_variables[i]) + (1|year_factor), family = binomial, data = spp_master_ztemp[reg == "newf"])
    bottom_extinction_model_comparison_newf[i,variable := bottom_variables[i]]
    bottom_extinction_model_comparison_newf[i,coef := coef(summary(mod))[,"Estimate"][2]]
    bottom_extinction_model_comparison_newf[i,p_value := coef(summary(mod))[,"Pr(>|z|)"][2]]
    bottom_extinction_model_comparison_newf[i,AICc := AICc(mod)]
    
    print(paste(i, length(bottom_variables), sep = "/"))
      
  }

#add column for those that didn't converge
bottom_extinction_model_comparison_newf$converge <- T
bottom_extinction_model_comparison_newf$converge[c(1:124)] <- F

#add column to designate region 
  bottom_extinction_model_comparison_newf$reg <- "newf"


#shelf (none didn't converge)

bottom_extinction_model_comparison_shelf <- as.data.table(matrix(nrow = length(bottom_variables)))
bottom_extinction_model_comparison_shelf[, variable:=as.factor(V1)][, coef:=as.numeric(V1)][, p_value:=as.numeric(V1)][, AICc:=as.numeric(V1)]
bottom_extinction_model_comparison_shelf[, V1 := NULL]
  
  for (i in 1:length(bottom_variables)){
    mod <- glmer(now_ext ~ get(bottom_variables[i]) + (1|year_factor), family = binomial, data = spp_master_ztemp[reg == "shelf"])
    bottom_extinction_model_comparison_shelf[i,variable := bottom_variables[i]]
    bottom_extinction_model_comparison_shelf[i,coef := coef(summary(mod))[,"Estimate"][2]]
    bottom_extinction_model_comparison_shelf[i,p_value := coef(summary(mod))[,"Pr(>|z|)"][2]]
    bottom_extinction_model_comparison_shelf[i,AICc := AICc(mod)]
    
    print(paste(i, length(bottom_variables), sep = "/"))
      
  }

#add column for those that didn't converge
bottom_extinction_model_comparison_shelf$converge <- T

#add column to designate region 
  bottom_extinction_model_comparison_shelf$reg <- "shelf"


#sa (12, 31, 35, 55, 63, 72, 82, 87, 92, 96, 123 didn't converge)

bottom_extinction_model_comparison_sa <- as.data.table(matrix(nrow = length(bottom_variables)))
bottom_extinction_model_comparison_sa[, variable:=as.factor(V1)][, coef:=as.numeric(V1)][, p_value:=as.numeric(V1)][, AICc:=as.numeric(V1)]
bottom_extinction_model_comparison_sa[, V1 := NULL]
  
  for (i in 1:length(bottom_variables)){
    mod <- glmer(now_ext ~ get(bottom_variables[i]) + (1|year_factor), family = binomial, data = spp_master_ztemp[reg == "sa"])
    bottom_extinction_model_comparison_sa[i,variable := bottom_variables[i]]
    bottom_extinction_model_comparison_sa[i,coef := coef(summary(mod))[,"Estimate"][2]]
    bottom_extinction_model_comparison_sa[i,p_value := coef(summary(mod))[,"Pr(>|z|)"][2]]
    bottom_extinction_model_comparison_sa[i,AICc := AICc(mod)]
    
    print(paste(i, length(bottom_variables), sep = "/"))
      
  }

#add column for those that didn't converge
bottom_extinction_model_comparison_sa$converge <- T
bottom_extinction_model_comparison_sa$converge[c(12, 31, 35, 55, 63, 72, 82, 87, 92, 96, 123)] <- F

#add column to designate region 
  bottom_extinction_model_comparison_sa$reg <- "sa"


#wctri (1, 2, 3, 5:10, 15:20, 25:29, 36, 38:42, 44, 47, 49, 52, 54, 55, 56, 58:65, 68:75, 79: 83, 85:87, 89, 90, 94:98, 100, 102:106, 109:111, 113:116, 119, 121, 124 failed to converge)

bottom_extinction_model_comparison_wctri <- as.data.table(matrix(nrow = length(bottom_variables)))
bottom_extinction_model_comparison_wctri[, variable:=as.factor(V1)][, coef:=as.numeric(V1)][, p_value:=as.numeric(V1)][, AICc:=as.numeric(V1)]
bottom_extinction_model_comparison_wctri[, V1 := NULL]
  
  for (i in 1:length(bottom_variables)){
    mod <- glmer(now_ext ~ get(bottom_variables[i]) + (1|year_factor), family = binomial, data = spp_master_ztemp[reg == "wctri"])
    bottom_extinction_model_comparison_wctri[i,variable := bottom_variables[i]]
    bottom_extinction_model_comparison_wctri[i,coef := coef(summary(mod))[,"Estimate"][2]]
    bottom_extinction_model_comparison_wctri[i,p_value := coef(summary(mod))[,"Pr(>|z|)"][2]]
    bottom_extinction_model_comparison_wctri[i,AICc := AICc(mod)]
    
    print(paste(i, length(bottom_variables), sep = "/"))
      
  }

#add column for those that didn't converge
bottom_extinction_model_comparison_wctri$converge <- T
bottom_extinction_model_comparison_wctri$converge[c(1, 2, 3, 5:10, 15:20, 25:29, 36, 38:42, 44, 47, 49, 52, 54, 55, 56, 58:65, 68:75, 79: 83, 85:87, 89, 90, 94:98, 100, 102:106, 109:111, 113:116, 119, 121, 124)] <- F

#add column to designate region 
  bottom_extinction_model_comparison_wctri$reg <- "wctri"

Now, I’ll merge all bottom extinction models

bottom_extinction_model_comparison_byregion <- rbind(bottom_extinction_model_comparison_ai, bottom_extinction_model_comparison_ebs, bottom_extinction_model_comparison_gmex, bottom_extinction_model_comparison_goa, bottom_extinction_model_comparison_neus, bottom_extinction_model_comparison_newf, bottom_extinction_model_comparison_sa, bottom_extinction_model_comparison_shelf, bottom_extinction_model_comparison_wctri)

Best bottom extinction models per region